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Abstract 



The paper presents a numerical technique for computing directly the Takens-Bogdanov 
points in the nonlinear system of differential equations with one constant delay and two 
parameters. By representing the delay differential equations as abstract ordinary differen- 
tial equations in their phase spaces, the quadratic Takens-Bogdanov point is defined and a 
defining system for it is produced. Based on the descriptions for the eigenspace associated 
with the double zero eigenvalue, we reduce the defining system to a finite dimensional al- 
gebraic equation. The quadratic Takens-Bogdanov point, together with the corresponding 
values of parameters, is proved to be the regular solution of the reduced defining system 
and then can be approximated by the standard Newton iteration directly. 

Keywords:Takens-Bogdanov point; delay differential equations; defining system; New- 
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; 1 Introduction 

■ In this paper, we discuss the computation of the Takens-Bogdanov points in the following delay 

| differential equations (DDEs) 

x(t) = f(x(t),x(t-r),X,n), (1) 

^ ■ where A, fi € R are parameters, r > is a constant delay. The DDEs have been widely studied 

because they often give more accurate descriptions for the phenomena in nature and engineering 
by taking into account not only the present state but also their histories (see [I], [7J, [2], [2] 
etc.). 

Takens-Bogdanov (short for T-B) bifurcation is one of the important bifurcations in dy- 
namical systems, it explains the mechanics for the occurrence of the Hopf point branch and the 
homoclinic branch as well as the saddle-node point branch. It acts as a bridge for investigating 
the global dynamical behavior through the local properties of the dynamical systems. Theoret- 
ically speaking, based on the local bifurcations near T-B points, one can obtain the T-B points 
on the paths of fold points (when a second real eigenvalue crosses the imaginary axis) or on 
the paths of Hopf points (when the two conjugate pure imaginary eigenvalues coalesce). For 
details, one can refer to [IS], [5] for ODEs and to [5], [H] for DDEs. 

The paths of fold points or Hopf points are generally traced numerically by the continua- 
tion method pQ. However, when T-B point is encountered along the paths, the continuation 
technique is not valid any more. As a result, to compute directly the T-B point is of great 
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importance. In addition, it is the foundation of branch switching, that is to compute the other 
solution branch, i.e., homoclinic branch or another solution branch differing from the paths 
traced before, emanating from T-B point. In general, the T-B point is approximated by numer- 
ical iteration method applied to an enlarged system. In this process, an initial value is necessary 
and can be obtained by detecting the changes of eigenvalues along the paths of fold points or 
Hopf points by the continuation technique [12]. For DDEs, the work could be done by the 
method developed by Luzyanina and Roose |16j , where a defining system is given to determine 
the Hopf bifurcation point algebraically, and the continuation techniques are considered. These 
techniques are enclosed in DDE-BIFTOOL [8]. Of course, to start the continuation of the path 
of fixed points, a good approximation is needed. In fact, in terms of that the constant delays do 
not affect the position of equilibria of the system, the continuation methods for the equilibria 
of ODEs could be employed to DDEs [T6] [17] directly. Unfortunately, noting that T-B point is 
a singularity of codimension 2, the techniques for computing the Hopf points developed in |16] 
for DDEs could not be applied to compute the T-B points directly. Besides, many applications 
require us to obtain the T-B points numerically first so that we can start our continuation from 
the T-B point (e.g. |10jV In addition, tracing the T-B points branch emanating from the singu- 
lar point of higher order needs to compute T-B points first as well. Therefore, to determine the 
T-B points numerically is an important task in the numerical analysis for dynamical systems. 

It is well studied for ODEs to calculate the T-B points, one can refer to Griewank and 
Reddien [13], Govaerts [11] . Beyn [5], Yang [20] and the subsequent articles by many authors. 
The method is described sketchily as follows. 

Consider the following parameterized dynamical systems 

x = f(x,\,n), ieM", (2) 

where A, fx S M are parameters, / : M™ x M x R — > R™ is a continuously differentiable function 
with f(x°, A , = 0, i.e. x° is an equilibrium of ([2]) for A — A and \i = 

Definition 1 20 (x°, A°,/i°) is called a T-B point of {l]) if the following conditions hold: 

(i) /(x°,AV°)=0; 

(ii) Af(f%) = span{r? }, m + 0; 

(iii) K(f>) = {,eW\$, = 0}; 

(iv) %eft(/°). 

The conditions in Definition [T] imply that there exist r/i S M and £i £ Mi , such that 

fSvi + Vo = 0, 
lo Vi = 0, 



and 



(/ x °) r a+£o = o, 

£^0 = 0, 



where R n = 7V(/°) © M, M™ = A/"((/°) T ) © Mi, ^771 = is equivalent to m £ M, ([uj q = is 
equivalent to £1 S Mi. 



Definition 2 [20] (x°, A°,/i°) is called a quadratic T-B point of (0), if it is a T-B point and 
fx 



(i) £o T / A ° + 0; 



(ii) d = det _ 0; 
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(iii) ± o, 




The defining system for the quadratic T-B points of ^ is then given by 



H 1 (w 1 )= f x {x,\,u)Q + rj = 0, 



(3) 



V 'o T c / 



where Wi = (x, rj, (, A, n) T £U = W l x K n x 1" x E x 1, ffi : [/ -> [/. 

Theorem 1 [20] Assume (x ,A°,/i ) is a quadratic T-B point of t/ien Hi(wi) is regular 



By Theorem [Q the quadratic T-B point of © can be obtained numerically by applying 
iteration methods to the defining system <j3j) . However, unlike ([3]), the defining system for T-B 
point of (JTJ) should be defined in a Banach space from the viewpoint of [13], see section 2 for 
details. If one wants to solve such a system directly, certain discretization must be applied and 
the discretization error cannot be avoided. 

In this paper, we present a numerical technique for computing the T-B points in nonlin- 
ear system of DDEs with one constant delay and two parameters, which carries the results of 
computing the T-B points of ODEs into the case of DDEs. Based on the descriptions for the 
eigenspace associated with the double zero eigenvalue, we reduce the defining system to a non- 
linear algebraic equation with finite dimension to avoid discretizing it directly. The quadratic 
T-B point, together with certain values of the parameters, is proved to be the regular solution 
of the reduced defining system. Therefore, the quadratic T-B point of DDEs can be solved 
numerically by the classical iteration methods. 

The paper is arranged as follows. Through representing the DDE as an abstract ODE in 
Section 2, we define the quadratic T-B point and produce a defining system for the quadratic 
T-B point of DDEs according to the methods described above. In Section 3, by using of the 
descriptions for the eigenspace associated with double zero eigenvalue, we simplify the definition 
of the quadratic T-B point of DDEs, meanwhile, the defining system obtained in Section 2 is 
reduced to a finite dimensional one, which is proved to be regular at the quadratic T-B points 
such that can be solved by the standard Newton iteration. To show the efficiency of our method, 
a numerical experiment is carried for a predator prey system with time delay in the last section. 

2 Regular defining system for the T-B points in Banach 
space 

By a simple change of timescale, one may take r = 1 in ([!]). Hence we consider the following 



where X,fi £ K are parameters, x £ R™, f(x, y, A, /i) is a C r (r > 2) smooth function from 
1" x 1™ x M x I to M". Assume x° is the equilibrium of f(x(t),x(t — 1), A, fi) at parameter 
value A = A and u, = u°. Eq. (j4]) can be linearized at (x°, A , as 



at u>i 



,9 - 



(x°,n , m ,\°,v°) T . 



DDEs 



x(t) = f(x(t),x(t-l),\,n), 



(4) 



x(t) = f°(x(t) 



x°) + f$(x(t-l) 




(5) 
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where /° = ^(x°,x°, A , ^°), ft = % (x°, x°, A , 

Denote by C := C([— 1, 0], M") the Banach space of continuous mappings from [—1, 0] to R™ 
with norm = maxe g [_ 10 ] |<H$)| (| • | is a some norm in R n ), and let 

^(9) = { f 2 (x ,x°,\,n), #e (-1,0), 

I 0, 9 = -l 

be a bounded variation matrix-valued function on [—1,0] parameterized by A and [i. Noting 
that 

o 

dr) X>fl (0)(x t (6) - x°) = fi(x°,x°, A, u.)(x(t) - x°) + f 2 (x°, x°, A, n)(x(t - 1) - x°), 

-l 

where x t (Q) — x(t + 9), the linear operator from C to R ra defined by 

L\^x t = / dn\^(9)x t (9) 



is bounded. 

The solutions of §5§ generate a Co-semigroup {T (t), t ^ 0} on G with infinitesimal generator 
A :C — > C defined by 

Ao4> = 0, 

D(A)= {^C 1 ([-l,0],r);^(0) = y ^,,.(^(0)}. 

It is known that the spectrum of the operator Ao consists of its point spectrum only [14j . i.e. 
cr(Ao) = <jp(Ao), and z S <t p (.4o) if and only if 

det(zl - /° - / 2 °e- z ) = 0. (6) 

Usually we call © the characteristic equation of ([S]). 

To develop the theory for computing the T-B points of DDEs, it is convenient to write ^ 
as the following abstract ODE in C [14] 

—u = G(u,X,fi), (7) 

where 

ru, \ n\(ff\-l /(«(o),u(-i),a,aO, ^ = o, ^ 

G ^ X '^ d >-{ u(9), 9 €[-1,0] (8) 

and the domain of G(-, A, /i) is {u 6 C 1 : u(0) — f(u(0), u{— 1), A, /i)}. 

Let C* = C([0, l],R n *) be the adjoint space of G, with R n * the n-dimensional space of row 
vectors. The adjoint bilinear form on C* x C is defined by [H] 

ty, 0) = V(0)^(0) - J 1>(£ - 0)d»/ A o i|t o (fl)^(Ode- 

The T-B point of ^ can be determined through investigating the T-B point of since they 
share the same equilibria. Here we define the T-B point of ([7]) first. 

Definition 3 (u , A ,fJ, ) is called a T-B point of ffl, or equivalently of Q), if 

(i) G(u°,A ,/i°)=0; 

(ii) Ar(GO) = span{0 1 },0 1 ^O; 

(iii) = {v e G, (^,«) = 0}; 

(iv) 0! e^(G°J. 
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The condition (iv) in Definition [3] is equivalent to (tp2,<Pi) = 0, and there exist ipi G G* and 
4>2 £ C such that 

(G°)*Vi+^2 = 0, (^i,wd) = 0, 
G°0a + & = O, (*o,0 2 ) = O, 

with wo G C and lo S G* satisfying 

(^2,W )-1 = 0, (f ,^l)-l = 0. 



(ii) do = dct 



7^0; 



Definition 4 (u , X ,fjr) is called a quadratic T-B point of or equivalently of ffl), if it is 
a T-B point and 
(i) W> 2 ,G° A )^0; 

(ip 2 , A4> 2 ) + (ipi,Act>i) (ip 2 , B4> 2 ) + (ipi,B(l>i) 

(iii)(^2,<fe) ^0, 

where A = G° uu (pi, B = G a au v Xll + c Aai G° a + G° M ura'i/i c Am = - (^'g») and Vx^ satisfying 

G u u\n + c\^G\ + Gj^ = 0, 
{lo,v\n) = o. 

Based on the preparations above, we produce the following defining system for the quadratic 
T-B points of gj) 

/ G(u,X,fx) \ 
G„(u, A,/i)ei 

H 2 {w 2 )= G u (u,\,n)e 2 + ei = 0, (9) 
(io,ei)-l 
V (io,e a ) / 
where w 2 = (u, e x , e 2 , A, e F = G x G x C x K x E, H 2 ■ Y -> F. 

Theorem 2 Assume [u ,X ,fi ) is a quadratic T-B point of Q), or equivalently of |^), i/ien 
^2(^2) is regular at w 2 = (u°, <f>\, (f>2, A , n°) T ■ 

Proof. The proof is quite similar to Theorem [T] □ 



3 Simplification of the defining system 

The defining system ^ is defined in a Banach space, so if we solve it directly, we will encounter 
many difficulties. The space G must be discretized first whenever what kinds of numerical 
methods are applied. This needs to store a large amount of data and causes the discretization 
error undoubtedly, which is not expected surely. Another big difficulty lies in the form of the 
function G(u, A, /i) which needs to be dealt with piecewisely since its special form caused by 
representing the DDE as an abstract ODE. 

In this section, we will reduce © to an equivalent form, which has finite dimension and can 
be solved easily. 

We first give an equivalent definition for the T-B points of (0} by the infinitesimal generator 
Ao of the Go-semigroup defined by the solutions of ([5]). 

Definition 5 A°,^°) is called a T-B point of Q), if 

(i) f(x°,x°,X°,^)=0; 

(ii) Af(Ao) = span{0i},0i ^ 0; 

(iii) K(A Q ) = {ieC, (V> 2)7 ) = 0}; 

(iv) 0! eK(Ao). 
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i) = 0, and there exist (p 2 G M and 
= 0, 



G M with I = h-sl 2 , 



The condition (iv) in Definition [5] is equivalent to (tp 2 , 1 
tpi G Mi such that 

A02 + 01=O, (Z,0 2 ) 
-4^1+^2=0, (-01,TU)=O, 

where C = Af(A ) ®M, C = Af(Ao) © Mi, (Z, <fe) = is equivalent to 
s G [0, 1]; (^i,ro) = is equivalent to ipi G Mi. 

Noting that x° G R™ can be regarded as a special element in C coincided with u°, (x°, A , uP) 
is a quadratic T-B point if it satisfies Definition [5] and |H The following theorem can be used 
to determine whether (x° , X° , n°) is a T-B point of ^ or not. 

Theorem 3 [19j Assume Rez ^ if z E <j p (Ao)\{0}, Eq. ftas a T-f? singularity if and 
only if the following conditions hold: 

(t) rank(/ 1 ° + / 2 °)=n-l; 

(ii) i/AA(A° + /°) = span{0?}, tten (/° + 7)0? G + / 2 °); 

(m) */ (a° + / 2 °)$ = (/ 2 ° + (/ 2 ° + ^ - \m t nn + / 2 °), 



where 



Denoted by P the invariant space of A® associated with the eigenvalue zero and P* the dual 
space of P, $= (<f>i(9),fa(6)),-l < < and V = col(ipi(s),i> 2 (s)),0 < s < I the bases of P 
and P* correspondingly, we have ($',<£) = I. Moreover, the following lemma holds. 

Lemma 4 |19] The bases of P and its dual space P* have the following representations: 

P = span<Z>, ${9) = {<f> 1 (6),<f> 2 {9)), -1 < 9 < 0, 
P*=spantf', gq>) =col(ipt(s),ip2(s)), < s < 1, 

w/iere <^(0) = 0° G R" \ {0}, ^C^) = 0§ + $9,$ G K" and ip 2 (s) = ip 2 G M' 1 * \ {0},^i(s) = 
■0? - sip 2 , V'i G K"*, w/iic/i satisfy 

(1) (A° + iM = o, 

(2) (/i° + / 2 O )^ = (/ 2 °+/)0?, 

(3) ^ 2 (/i + / 2 °) = 0, 

(4) <(A + / 2 ) = ^ (/ 2 ° + /), 

(5) ^i-u° 2 m+^m = h 

(6) - him + + - him = o, 



(10) 



where we can determine the unique vector <p{ , ip^ by (1) and (3), respectively, up to some constant 
factors; then we can determine 4> 2 ,tpi by (2) and (4), respectively. However (5) and (6) are 
used to determine the coefficient factors of the vectors (fr® and ^ 2 . 



Based on the characterizations above, we can give another definition for the quadratic T-B 
points of ([4}. 

Definition 6 {x , A ,/x ) is called a quadratic T-B point of ^jj, if it is a T-B point and 

^(B 1 + B 2) 

V> 2 °(P1 + P 2 ) 



( ^{A x + A 2 )0? 



(ii) do = dct 



(iii) V 2 >§ 



where 



^ 1 {A 1 +A 2 )<f>\+ 
^ 2 {A X +A 2 )^- 



^0; 



^i = (/i°i+/i° 2 ) 

Bi = (A°i + /i>A M + cxpfa + ftp B 2 = {f 2 \ + / 2 °> v + CAM / 2 ° A + / 2 % 



^2 = (/°1 + / 2 ° ; 
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with c Xtl — -^rjv and v Xll satisfying 

(A° + f>x, + cx m /a° + ft = 0, 
^x >fl = 0. 

Theorem 5 Definition^ and\^ are equivalent. 

Proof. We only need to show that the corresponding conditions in Definition [4] and [6] are 
equivalent. We show the conditions (i) in both Definition |4] and [6] are equivalent first, that is, 
(ip 2 , G°) ? V 2 °/a ± 0. Noting that 

G ° = { o*' leKO) 

and 4 x -u{6) — as e [—1,0), we obtain 

(V2,g a ) = v 2 (o)/^- /° [ M^-o)[dvxo,Ae)]G x (u ,\°,fjP)(OdC = ^fl (12) 

7-1 Jo 

Besides, ("02, ^2) = ^202 — l^ss/aV? + ^2/2 $2 > this confirms that the conditions (hi) in both 
Definition |4] and |6] are equivalent. We only need to show that the conditions (ii) in both 
Definition |4] and |6] are equivalent next, and it will be a tedious calculation. Firstly, we get 

= ^ 2 °(/iVi(0)<M0) + A 2 ^(o)^(-i) 

= V> 2 °(AVi>? + A° 2 « + / 2 °i^S0? + A° 2 «) ( ' 

= M& + 

In a similar way, we have 

= + A°2 + /2°1 + 

+^2°(/i°A + / 2 A)CA^? + V 2 °(A°„ + n 4] 

= ^[(/l l + /i°2+/2 l+/2 2K, 

+ (A°A + /2 A)^ + (A /i + /2° M )]^ 

Since f® 2 = / 2 i, we obtain 

(V< 2 ,^ 2 ) =(^ 2 ,G° U ^^) 

= (l/&(GM)($ + #0)) 

- i$[/iV2($ + 0? • 0) + aV?($ + $ • -1) 
+/ 2 Vi(02* + 0? ■ 0) + A° 2 0i(^ + 0? • -i)] 

= + A° 2 0?($ - + / 2 ViV 2 + A° 2 0i(^ - 0?)] 

= ^2°(A°i + A° 2 + A°i + f^m - ^ 2 °(A°2 + 

=v 2 °(^i+a 2 )^-v 2 °^;. 

Jointing the above equation with 



= (^-^ 2 °,G° u «) 
= + A 2 )</>°, 
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we obtain 

(V>2, Afa) + (^i, A<fn) = V 2 Vi + A 2 )(b° 2 - tfiAtfl + ^(A! + A 2 )cb° 1 . 
At last, noting that 

{ih, Bfo) = (V2, (G^ + c x ,G° uX + Gl^ 2 ) 

= m>xM + 4>l ■ 0) + f>xM + 4>\ ■ -i) 
+f>xM + 0? • 0) + / 2 ° 2 ^(^ + 0? • -i) 
+ca„/ia($ + 0? • 0) + c v /° A ($ + 0? • -1) 
+/i O ^(^ + ^-O) + / 2 ° Al (^ + 0?--l)] 

= V^I/iVa^ + f&\M - 0?) + / 2 >A^ + / 2 ° 2 ($ " 0?) 

+cam/?a^ + + M + f§M - ft)] 

= ^ 2 {B l+ B 2 )4>l-4> Q 2 B 2 ^ 



(15) 



and 



<[(Ai + A°2 + /: 

f°A 



+^(/ 1 ° A +/ 2 O A )0? + (/ 1 O A1 + / 2 ° AI )0?] 
AO 



we have 



(V>2 



2 ) + {faBfa) = $(Bi + B 2 )cj>° 2 - ^B 2 <p° + i>l(B t + B 2 



(16) 



Collecting Eqs. (|T3|) to (fT6|) . we arrive at our assertion. □ 

According to Theorem [3] and Lemma |4l we introduce the following defining system for the 
quadratic T-B points of (@| 



H(v) 



( f(x,x,X,n), 

(fi(x, x, A, n) + f 2 (x, x, A, n))tpx 

(h(x, x, A, n) + f 2 (x, x, A, n))ip 2 - (f 2 (x, x, A, n) + I)(pi 



V 



hfi - \hh{x, x, A, u)<pi + hf 2 (x, x, A, (i)<pi - 1 
Ziy> 2 - 3/1/2(1, x, A, n)tpi + hf 2 (x, x, A, n)<p 2 + 



= 0, 



(17) 



ll 2 f 2 (x, x, A, jti)^i - 5/2/2(2:, cc, A, |«)</J2 / 



where v = (x, (p ll <p 2 , A, fi) T e V = M" x R n x 1 
above, we know the system is equivalent to Q. 



°, H : V — > V. From the discussions 



Theorem 6 Assume (x° , A , u°) is a quadratic T-B point of then the defining system |7 
is regular at its zero v° = (a; , <f>\, <p 2 , A , fi Q ) T ■ 
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Proof. Obviously, H° reads 

/ /i° + / 2 ° fO f° 



-^i(/i° 2 +/ 2 ° 2 )0?+ ri^ ?+ 

V I'2(/i° 2 +/ 22 )^ \hf^l 

We first prove the map H® is injective. By expanding H®§ = with § = ($1, $ 2 , ^3, ci, C2) T G 
V, we obtain 

(/i° + / 2 >i+ C i/ A ° + C2/; = 0, (18) 
(A°i + A° 2 + / 2 °i + &)<t>>i + (A° + / 2 >2 + {fix + / 2 °a)^?ci + {fl + fljfici = 0, (19) 
\fn + /ft + / 2 °i + - (/a ! + / 2 ° 2 )<£?)tfi - (/ 2 ° + T)* 2 + (A + / 2 °)tf 3 

+((a°a + / 2 °a)^ - /m)* + ((/?, + - fUih = °- 



( 



and 



5fe(/£> + / 22 )^i + MA° 2 + /Mj* + (*i - \hft + hf$)#2+ 

+ (-3fe/lV? + + i-^fUl + 'l/2>I>2 = 

M? a +/M+ll(/ft+/s^° 



(20) 
(21) 



( 



2 1 H./12 ~ J 22 m i >HJ12 -r j 22 /v 2 

+|MA° 2 + - §fc(/? 2 + f§2)$)4l 

+(-5*1/2° + \hf^2 + (h + hfl - |/ 2 / 2 °)^3 (22) 

+ (-3'l/ 2 A^l + l ^ 2X <f 2 + |*2/ 2 °A - Ihf^C! 
H-^hfUl + l lfU°2 + Sfe/2° M ~ 5*2^)02 = 0. 

To multiply (fT"8| by ip 2 from left, we have by Lemma |4] 

^/aCi+^/°c 2 = 0, 

which yields c\ = — ^pjt>c 2 . Noting that c\^ = — ^„^S , we obtain c\ = c\^c 2 , therefore, 

$i = c(f>1 + c 2 v\^, where c is a constant to be determined. By substituting them into Eqs. (|19[) 
and (fJDD, it yields 

c(Ai + A 2 )<jP 1 + c 2 (B 1 + 5 2 )0? + (A + / 2 °)tf 2 = 0, (23) 

and 

c(A 1 + A 2 )^ + c 2 (B 1 + B 2 )<Pl - cA 2 $ - c 2 B 2 <\>\ - (/ 2 ° + I)-d 2 + (f° + / 2 > 3 = 0. (24) 
To multiply (1231) by ip 2 from left, we obtain by Lemma 0] again 

ct/^i + A 2 )0° + c 2 V 2 (Bi + B 2 )0i = 0. (25) 
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To multiply ((23]) by ipi from left, we have 

oiPliAx + A 2 )<j>\ + c 2 ^(Bi + S 2 )0? + + / 2 >2 = 0. (26) 

Multiplying ((24]) by from left yields 



-c^\B^\ - ^ 2 °(/ 2 ° + I)d 2 + + / 2 °)i? 3 = 0. 

Adding (f26| and (f27|) together, and utilizing (|10|) we obtain 

c(^ 1 (A 1 + A 2 )^ + </>°(Ai + A 2 )$ - $A 2 $) 

+c 2 (^?(Bi + B 2 )<^ + ^{Bx + B 2 )<p° 2 - ^Bxft) = 0. 



(27) 



(28) 



From ((25)) and ((28]) , noting that do ^ 0, we obtain c = c 2 = 0, therefore c\ = Q, #i = 0. 
Consequently, ((Ml) reads (/£ + / 2 °)tf 2 = 0, therefore, ?? 2 6 A^/? + J 2 °). If i? 2 7^ 0, we have 

h$2 ~ \h!l^2 + h.f^2 = M2) ± 0, 

which contradicts with (|2~1T) . hence i9 2 = 0. Hereby, ((2"0"|) reads (/° + / 2 )$3 = 0, together with 

d22) it implies tf 3 = 0. Altogether d = 0. 

In a similar way, we can prove H® is surjective. This concludes the theorem. □ 

Noting that the system ((17]) is of finite dimension, according to Theorem (6( it can be solved 

by many iteration procedures, especially the standard Newton iteration procedure, that is 

v k +i = v k - [Jnivk^Hivk), for k ^ (29) 
where Jh(') is the Jacobian matrix. 

4 Numerical example 

Consider the following predator prey system with delay |15j . 

\ i 2 (t) = X2 (t){J^0L - D) [ 1 

where r, K, a, fx, D and r are positive constants. Eq. (|30[) has a T-B singularity if the parameters 
/1, a, D and if satisfy u 2 — 4aD 2 = and /1 = ifD |15j . For determining the T-B point of (I3"0"l) . 
we choose if and D as parameters and fix the others. Specially, taking o = 1, /z = 1, r = l, 
and r = 1 in (|30p . we have 

i(t) = /(a;(t),a:(t-l),A*') (31) 

with x(i) = (xi(t), x 2 (t)) T and 

/ a:i(*)(l - ^#) - ^".V* 2 ^ \ 

V z a(*)( 1+ ^ (t _' l 1) - j>) ; 

For any given iel 2 , the derivatives of / with respect to the first and the second variables at 
x are respectively given by 



h(x,x,D,K) = 



-1 2xi x\ 

1 K T+xJ 

-rp-T -D 

l+xf 
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and 



f 2 (x,x,D,K) 



( Tn £=gL n 

\ X2 JT+^f U 



Hence the denning system for the quadratic T-B point of ([50)1 reads 
/ f(x,x,D,K) 

{f 1 {x,x,D,K) + f 2 (x,x,D,K)) l p l 
(AO, x, D, X) + / 2 (x, D, K))tp 2 - (f 2 {x, x, D, K) + 
hipi - \l 2 f 2 (x,x,D,K)(pi + hf 2 (x,x,D,K)ipi - 1 
h<p 2 - \hf 2 {x, x, D, K)ipx + hf 2 (x,x, D, K)ip 2 + 



H(v) 



\ 



0, 



(32) 



hhf 2 (x, x, D, K)ipi - \l 2 f 2 (x, x, D, K)ip 2 I 



where v = (x, <fi,ip 2 , £>, K) T E V = K 2 x K 2 x R 2 x M x E, i? : V" -)• K. 

Choosing Zi = (1,0) and / 2 = (1, 0) and applying Newton method (|2U)) to Eq. (f3"2"j). we 
obtain the following results shown in Table 1 with Matlab. In each iteration with respect to 
different initial values, the solution converges to (1, 1, 1, 0, 0, —2, 0.5, 2) with the remainder less 
than 10 -21 . 



Initial value vo 


Step 


(1.1, 1.1, 1, 0, 3, 0, 0.4, 1) 


5 


(1.2, 1.2, 1.2, 1, 1, 0, 0.5, 0.5) 


7 


(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 0.6, 1.6) 


7 


(3, 1.5, 1.2, 0.5, 1.8, -1.8, 0.45, 1.9) 


6 



Table 1: Number of iterations required by Newton methods for (|32|l with respect to different initial 

values. 

In fact, according to [T3] we can obtain exactly the T-B point of (|3U1) . For o = /U = r = r=l, 
the T-B point is (xi,x 2 ) = (1, 1) with parameter values D — i and K = 2. 

Remark 4.1 From the numerical result we see that the Newton's method converges very rapidly 
once the initial value is close to the true solution. However, choosing the initial point Vq to 
ensure a convergent Newton iteration for solving J_?7| j is a rather complicated problem because 
the function H(v) is highly nonlinear and high dimensional. There is no good method as so far 
even for ODEs or for the first order singular points without the aid of continuation techniques. 
The guess of the initial value Vq for starting the iteration procedure can be found in the following 
way. First we fix one parameter, for instance fi — u. H . According to the method by Luzyanina 
and Roose we can find a Hopf bifurcation point x H of on its Hopf bifurcation curve 
at A = X H , as well the approximation of the eigenvector (f>^ e + i<f>y m associated with eigenvalue 
ioj H . Next, we start from (x H , X H , [i H , lj h , 4>R e , 4>¥m) t° t™ce the Hopf bifurcation curves by 
the continuation technique. When the singularity of the nonlinear system used for continuation 
changed at some point, for instance (xq, Xq, [1q, ljq, <fi a Re , <ft] m ), then we can use (xq, <fi a Re , (j>, Xo, 
as initial value to start the Newton iteration procedure, where <f> still should be guessed, but it 
is not a hard work any more. 
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